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Abstract 

The past few years has brought renewed focus on the physics behind the class of materials characterized by long-range interactions 
*—* and wide regions of low electron density, sparse matter. There is now much work on developing the appropriate algorithms and 
,__! codes able to correctly describe this class of materials within a parameter-free quantum physical description. In particular, van der 
C^> Waals (vdW) forces play a major role in building up material cohesion in sparse matter. This work presents an application to the 
OQ vanadium pentoxide (V2O5) bulk structure of two versions of the vdW-DF method, a first-principles procedure for the inclusion of 
1— h vdW interactions in the context of density functional theory (DFT). In addition to showing improvement compared to traditional 
jj~ semilocal calculations of DFT, we discuss the choice of various exchange functionals and point out issues that may arise when 

treating systems with large amounts of vacuum. 
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1. Introduction 

Density functional theory (DFT) is one of the most reliable 
and widespread methods for atomic-level computational stud- 
ies of the structure and electronic properties of materials. The 
formalism behind DFT is exact, but for practical use an approx- 
imation for the exchange-correlation part E xc of the total energy 
functional must be found. 

Dense crystalline materials show two main characteristics: 
they have a periodic structure and they are characterized by 
strong bonds. The periodic structure simplifies the description 
in codes using periodic boundary conditions, and the strong 
bonds are well described already by one of the simplest of 
the approximations for E xc , the local-density approximation 
(LDA), and also by the later, semilocal, generalized gradient 
approximation (GGA). For this reason dense crystalline mate- 
rials have been studied systematically and in great detail for 
a long time. On the other hand the extension to account for 
matter with regions of low electron density and a considerable 
component of long-ranged forces, like the van der Waals (vdW) 
forces, has proven to be more difficult. The vdW interaction is a 
quantum-mechanical phenomenon providing bonding by corre- 
lating instantaneous charge fluctuations. A successful method 
that is able to describe the vdW forces, the vdW-DF method, 
has been developed relatively recently |T]|2l[3]. 

In this paper we present a DFT study using the vdW-DF 
method for vanadium pentoxide (V2O5). The a-phase of V2O5 
has a layered structure, illustrated in FigurefT] The bonds within 
the layers are strong while the interlayer binding is long-ranged 
and weaker and for this reason the material is easy to cleave in 
the plane perpendicular to c. 



V2O5 has previously been studied with DFT by a num- 
ber of groups using GGA H H IS or by adding 
(semi-)empirical vdW-terms to GGA calculations [9]. Here 
our focus is on the vdW-DF method rather than on the mate- 
rial V2O5 itself. We therefore mainly discuss the issues aris- 
ing when using vdW-DF in a layered, extended system such as 
V2O5 bulk. 

The structure of this paper is as follows. In Section 2 we ex- 
plain the computational methods used, focusing in particular on 
the technicalities required both by the self-consistent GGA cal- 
culations and by the vdW-DF postprocess procedure that allows 
the inclusion of vdW forces. In Section 3 we show and discuss 
our results, and Section 4 contains a summary. 




Figure 1: Illustration of the V2O5 bulk structure. The unit cell is here repeated 
four times to illustrate the layered structure. The lattice parameters a, b and c 
are indicated. The black spheres are the vanadium atoms and the small spheres 
are the oxygen atoms |10|. 



'Corresponding author 
Email address: londeroSchalmers (Elisa Londero) 



Preprint submitted to journal for publication 



July 20, 2010 



2. Calculational methods 

Vanadium pentoxide has a layered structure with charge 
voids between the layers, and for this reason a DFT analysis 
requires use of the vdW-DF method. The calculation consists 
in moving the vanadium pentoxide layers apart in the stacking 
direction, accordion-like. We study the changes in total energy 
with changes in the values of c. The intralayer binding, with 
a mixed covalent and ionic nature, is a strong binding and this 
justifies our choice of keeping the interatomic distances fixed 
inside the layer itself when the layers are taken apart. 

The computational procedure is made up by two steps. First 
a self-consistent GGA calculation using the plane-wave DFT 
code dacapo ifTTI is carried out and subsequently the vdW-DF 
method is applied to calculate the nonlocal energies in a post- 
GGA procedure. This has been shown to not produce signif- 
icant discrepancies in the energy calculations, compared to an 
entirely self-consistent vdW-DF calculation [2|. 

Focusing on the details of the self-consistent GGA calcula- 
tions, the Brillouin zone is sampled using a 2x4x4 Monkhorst- 
Pack £-point set and both the plane-wave and density cutoffs 
are set to 500 eV. We check the convergence of our calculations 
with respect to these parameters. The orthorhombic unit cell 
contains two V2O5 formula units and it is periodically repeated 
in the three spatial directions. We use in-plane lattice constants 
a - 11.55 A and b - 3.58 A, which are the optimal values in 
GGA calculations. Ultrasoft pseudopotentials (USPP) are used. 
The fast Fourier transform (FFT) grid is chosen such as to have 
at most 0.12 A between nearest-neighbor grid points. For each 
FFT grid point the electron density n(r) is self-consistently cal- 
culated within GGA and subsequently a vdW-DF calculation is 
performed starting from n(r). 

We take isolated layers of V2O5 as the reference point of the 
energy. The reference calculations are carried out using a unit 
cell that has a length in the c direction that is 4 times the length 
of the bulk unit cell, with a single layer of V2O5 in the middle 
of it, surrounded by vacuum. 

For the inclusion of the vdW interactions in our calculations 
we use the vdW-DF scheme implemented in a non-periodic 
code. The correlation energy E c is divided into two parts, one 
part that is mostly local and approximated by the LDA correla- 
tion energy £'Jr DA , and one part that includes the most nonlocal 
interactions and is written [II: 



E% 



W 4//w„(r ¥ (r,rW). 
The correlation term is thus 



F LDA + Ef 



(1) 



(2) 



and it is added to the total energy from the self-consistent GGA 
calculation after removing the GGA correlation E^ GA . 

As mentioned above, the code used for this second step is 
not periodic in space so the spacial periodicity of the material is 
reproduced by adding a number of cells in the three directions 
of space (for the bulk calculation). More precisely, for each 
electron density grid point r within the central unit cell, all the 
other grid points r' to be used within the integral (fTh fall in 



two regions, limited by two radii. Inside the smallest radius all 
available grid points are used in the Ef calculation. In the other 
region — inside the larger radius but outside the smaller radius — 
only half of the grid points in any of the three spacial directions 
are used 1121 . The reference calculations (isolated layers) are 
carried out similarly, with the exception that no additional grid 
points are included in the stacking direction c. The two radii are 
chosen such that further inclusion of grid points only changes 
the energy contribution from Ef marginally. 

At this point another computational choice must be taken: the 
most appropriate exchange functional to be used. The vdW at- 
traction is purely a correlation effect so binding from exchange 
that mimics the vdW binding must be avoided. The revPBE ex- 
change functional has often been used because for some exam- 
ples of sparse matter it has been shown to give the least spurious 
exchange contribution E x to the binding energy [1|. However, 
revPBE exchange is overly repulsive, and it is relevant to also 
consider other choices of exchange functionals. The exchange 
forms used here are discussed in Section 3. 

There are two relevant versions of the vdW-DF method 
presently available [J The first 03 0], which is here called vdW- 
DF1, has shown to work well by improving binding energies 
and separation distances over the GGA results. It has been 
tested lfl5l on a variety of systems among those: atomic and 
molecular systems (such as dimers of benzene, benzene-like 
molecules, and polycyclic aromatic hydrocarbons (PAH), and 
polymer interactions), crystalline solids (like graphite, potas- 
sium intercalation in graphite), and adsorption (benzene and 
adenine on graphite, PAHs on graphite, adsorption of aromatic 
and conjugated compounds on M0S2). Despite these successes, 
vdW-DFl underestimates the hydrogen-bond strength and over- 
estimates equilibrium separations. 

The slightly different version of the vdW-DF method, called 
vdW-DF2 [3 1, has recently been developed. It has been con- 
structed focusing on obtaining more accurate energies for finite, 
relatively small molecules, but here we test it for the extended 
system vanadium pentoxide. It uses a different expression of the 
plasmon frequencies used in the evaluation of Ef, To reach a 
description within chemical accuracy for most of the molecules 
of the S22 set, the use of a different semilocal exchange func- 
tional is also suggested in Ref. J3). The exchange choice sug- 
gested in Ref. [3 1 to go with vdW-DF2 is the refitted version of 

PW86 USEE). 

In the present paper we combine the two versions of vdW- 
DF with a number of exchange choices (subscript x denotes the 
exchange part of the functional): the revPBE v [ 18 1 and PW86 A 
choices mentioned above, and the recently suggested PBEsol. v 
Ifl9l . optB88 x GO), and C09* GD- 

3. Results and discussion 

We have previously presented V2O5 bulk binding results ll22l 
found by use of vdW-DFl. The choice of the exchange func- 



1 A preceding version H3II14I was developed for layered material in which 
the layers are approximately translationally invariant. This is not the case for 
V2O5 and thus that version is ignored here. 



tional, even if not the major step of the whole procedure, is of 
importance and the selection or the construction of the most 
appropriate exchange form is still a matter of debate. 

In this spirit we here use both vdW-DFl and vdW-DF2 in 
combination with five different forms of exchange, some of 
which for physical reasons may be considered good candidates. 
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Figure 2: Top panel: vdW-DFl energy curves for different choices of the ex- 
change functional compared with the results for vdW-DF2 (bottom panel). The 
vertical black lines show the optimal binding distance (no binding energy avail- 
able) according to experiments |23 1. The stars show the position of the mini- 
mum according to DFT calculations (GGA) that do not account for vdW inter- 
actions. 



3.1. Interlayer binding results 

Figureplshows the resulting total energy curves, £' vdw " DF (c), 
obtained from bulk calculations with the separated oxide lay- 
ers as reference calculations. As described in detail elsewhere 
1 22, 24 1 we ensure that the atomic positions on the underly- 
ing FFT grid in the reference calculation are identical to those 
of the bulk calculation. As discussed further below, we have 
excluded the contribution to Ef from a few grid points in the 
vacuum region that have unphysical values of n. The minima 
of the £ ,vdw-DF (c) curves (the binding energies Ef, with positive 
values for binding systems) are summarized in Table [Tl along 
with experimental results and a GGA calculation, using PBE 
|251. without inclusion of vdW forces. 
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4.72 


0.86 
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4.46 


1.19 
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4.28 


1.10 




PBEsol x 


4.28 


1.48 




optB88. v 


4.37 


0.39 


vdW-DF2 


revPBE, 


4.87 


0.48 




PW86 V 


4.55 


0.75 




C09 v 


4.38 


0.62 




PBEsol A 


4.36 


0.99 




optB88. t 


4.47 


-0.08 



Comparison PBE 4.87 0.18 

Experiment 3 4.368 



a Ref. [23]. 



Table 1: Results for the equilibrium distances and the binding energies per 
unit cell. The results are shown for both vdW-DFl and vdW-DF2 with various 
exchange functionals. The boldface entries are the defaults for the vdW-DF 
version. 



The various exchange choices of GGA type that we consider 
here differ by the choice of enhancement factor F x (s) in the 
exchange energy density. The exchange energy may be written 



E x = j dr «(r) e^ A («(r)) F x (s(r)) 



(3) 



with e A r DA ~ n l/3 and where s = c. s |Vn|/« 4 ^ 3 with c s - 
l/(27r 2 ^ 3 3^ 3 ) is the reduced electron density gradient. We plot 
F x (s) for various exchange forms in the insert of Figure [3] 
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Figure 3: Error measure rj for the contribution of pseudopotential-related 
"noise" to the exchange energy if not explicitly excluded. Insert: Enhancement 
factor F x (s) for the exchange energy. 



Several physical constraints may be imposed on the form of 
F x (s), as described in the literature [26|. In the limit of small s 
(approximately < s < 1), the functionals PBEsol v and C09 A 
are designed to follow the simple gradient expansion approx- 
imation (GEA) form. The small-i region corresponds to the 
region of large densities n. This was suggested in Ref. Ifl9ll to 
remove an artificial bias towards free atoms and in Ref. 12T1 



to reduce the short-range repulsion compared to revPBE A . The 
other three exchange choices (revPBE A , PW86 A , and optB88 A ) 
all have steeper growth of F x (s) for small s (Figure[3]l. 

In the large- s limit it has been argued ll27l IT7I that constant 
asymptotes of F x (s) may lead to the same type of spurious bind- 
ing from exchange alone as in LDA. This arises in regions that 
have both low values of the density n and an inhomogeneous 
distribution (large Vn), for example outside surfaces and in in- 
ternal voids. In Ref. ITTI it is argued that for small molecules 



a growth of - 



,2/5 



gives the best agreement with Hartree-Fock 



results. However, in our present system with extended oxide 
layers this is not necessarily the case. Of the exchange forms 
used here the revPBE A , C09 A , and PBEsol A tend to constants 
at large s, PW86 A grows as ~ s 2 ^ 5 , and optB88 A has a much 
steeper growth (~ s). 

Figure [2] and Table [T] show that in general, the improvement 
of the equilibrium distances using the vdW-DF technique is 
evident when compared to a standard GGA calculation. It is 
also clear that the suppression of F x (s) for small s in PBEsol A 
and C09 A results in smaller binding distances. Of all exchange 
choices presented here, PBEsol A and C09 A give the smallest 
binding separations. For vdW-DF2 those separations are very 
close to the experimental value, for vdW-DFl they are a few 
percent smaller than the experimental value. 

3.2. Numerical noise 

Soft matter have substantial regions of low electron density. 
In such systems the electron density can be as low as 10 -9 
|e|/A 3 . An issue related to the presence of such densities is 
that the reduced density gradient s ~ Vn/n 4 ^ that enters the 
enhancement factor F x (s) becomes extremely large unless the 
small density is locally almost constant (Vn small). Depending 
on the shape of F x (s), the exchange energy p) can be more or 
less influenced by those large s values. 

Our use of USPP and, more generally, noise in planewave 
code calculations, gives rise to an additional problem: the pres- 
ence of small unphysical negative electron density values. The 
dacapo code formally tries to solves this problem by introduc- 
ing a floor for the electron density, «fl 00 r, dac = 10~ 9 |e|/A 3 , and 
setting all the negative values equal to the floor before evaluat- 
ing the E x term. This procedure of course increases the number 
of very small but positive n values that give rise to very large 
values of s. 

In Figure |4] we show the distribution of s values for both a 
bulk and a reference calculation of the V2O5 structure. The 
plot is based on the valence electron density n from the dacapo 
GGA calculations. For evaluating s for the histogram in Figure 
|4]we have introduced the same type of replacement of negative 
n with a small positive value, since negative n make no sense in 
the evaluation of s. We here use the value nfl oor = 1CT 15 |e|/A 3 , 
but the distribution is not sensitive to the precise value of the 
«fi 001 . In the bulk calculation hardly any negative values of n 
appear (~ 30 grid points have negative values out of a total of 
~ 10 5 grid points), and in the reference calculation all values in 
the region-5 bin (s > 1212) originate in point that had n < «fl oor 
before the replacement. 
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Figure 4: Distribution of .s values for V2O5 at the equilibrium distance. Values 
are obtained from dacapo generated GGA based electron charge densities n. 
Top part: bulk. Bottom part: reference calculation (isolated V2O5 layer with 
much vacuum). The corresponding distribution for bulk of dense materials is 
usually in the range < s < 3. 



The small values of n, leading to relatively large values of 
s in dacapo E x calculations, are thus a mix of physically cor- 
rect small values and "noise" arising from the USPP descrip- 
tion of n. In order to limit contributions to E x from the low- 
density values (most of them unphysical), we carry out our 
own E x calculations, based on the dacapo valence density n 
but now with removal of all contributions from grid points with 
n < nfl 00r = 1CT 15 M/A 3 . To illustrate the effect of this removal 
we have calculated the quantity: 



'=^L*"^ DA( "">4^ 



(4) 



which is the difference (as expressed by A) between the integral 
of the exchange energy for the bulk structure and the integral for 
the reference structure. The r spans the density grid points that 
have a value of n below the floor nfl oor , i.e., the points we remove 
in our calculations. The error measure 77 includes a division by 
the volume of the (bulk) unit cell V - abc because more (phys- 
ical) vacuum area enters as the layers are moved apart, and this 
physical contribution is not to be seen as an error contribution. 
The results are shown in the main panel of Figure [3] 

We find that for each choice of exchange the measure 77 is 
approximately constant. This is in agreement with our expecta- 
tion that contributions to 77 come mainly from the vacuum part 
of the unit cell in the reference calculation, a region that grows 
approximately linearly with unit cell volume V and thus c. We 
find that for most of the exchange functionals the error measure 
is small, at approximately 0.4 to 1.1 meV/A 3 , but for optB88 A 
it is much larger, at 3.4 meV/A 3 . This means that unless we 
explicitly exclude those points, the optB88 A will contribute to 
the E x an erroneous 0.6 eV per unit cell at the binding distance. 

Noise on top of larger electron densities still affects our cal- 
culations. This is so even if we have, by the exclusion, removed 
an important part of the erroneous contribution. The error mea- 
sure should therefore only be taken as an indication of the sen- 



sitivity of the exchange choice to possible numerical noise in 
n. It is clear that optB88 A is much more sensitive that the other 
choices, to the extent that we will not include optB88 A in our 
work on vdW-DF development. Nevertheless, it is interesting 
to analyse the behavior of optB88 A to learn about the numer- 
ical challenges which exist for accurate calculations with the 
vdW-DF method. 

Why is the optB88 A so much more sensitive to noise than 
the other exchange choices we consider here? It is clear from 
the form of optB88 A that numerical noise in systems with large 
regions of low densities n will have large contributions to E x 
from these regions. The growth of optB88 A for large s is so 
steep, F x ~ s ~ n ^ (even if |Vn| is assumed constant at a 
finite value), that the factor ne^ DA ~ n 4 ^ 3 in the integrand will 
not suffice to damp the contribution. We have 



7 optB88 
J X, voids 



f drn(r)^ DA (n(r))FT B8S (s(r)) (5) 

Avoids 



Jvoi 



constant I dr]Vn(r)| 

/voids 



(6) 



and the latter expression is obviously very sensitive to noise. 
The F x for the other exchange choices grow less steep, or not at 
all, at large s, and the contribution from the void areas is thus 
damped by ne x DA . 

While the USPP is in our case likely the main origin of noise 
(of which we remove an important part explicitly), this is a 
problem not only for dacapo or other codes that use USPP. Most 
DFT codes were written at a time when small electron densities 
were not relevant and the codes are thus not optimized for han- 
dling regions of very low electron densities. 

We do not want to give recommendations of any particu- 
lar exchange form based on the one system (V2O5) studied 
here. However, we do show that calculations are sensitive to 
the choice. We also find that optB88 A is not a valid candidate in 
systems like the present both for physical reasons (F x (s) growth 
too steep as discussed in Ref. [ 17]) as well as for numerical rea- 
sons as discussed above. 

The effect of the noise on the binding energies found in Table 
[T] is severe for the case of optB88 A (estimated error 0.6 eV), 
but smaller for the other exchange choices (0.1 to 0.2 eV). It is 
important to note that these estimates are approximate, that they 
may include also physically correct small values of n (although 
this is not the case in this particular system), and that noise on 
the larger values of n is not taken into account. 

4. Summary 

The structural properties of a vanadium pentoxide bulk struc- 
ture are calculated within the vdW-DF method, testing both of 
the presently available version: vdW-DFl of Ref. [1 1 and vdW- 
DF2 of Ref. [ 3 1 . A number of forms of the exchange functional 
are tested. The interplanar distance, which is usually severely 
overestimated in traditional GGA calculations, is found to be 
closer to the experimental value using these methods. This con- 
firms the expectation that the interplanar bond in vanadium pen- 
toxide has a large vdW component. 
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